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ABSTRACT 



Context. A number of sytems of multiple super-Earths have recently been discovered. Although the observed period ratios are generally far 
from strict commensurability, the radio pulsar PSRB1257 + 12 exhibits two near equal-mass planets of ~ 4 M e close to being in a 3:2 mean 
motion resonance (MMR). 

I Aims. We investigate the evolution of a system of two super-Earths with masses < 4 M e embedded in a turbulent protoplanetary disk. The aim 
O ■ is to examine whether or not resonant trapping can occur and be maintained in presence of turbulence and how this depends on the amplitude 
l— > ' of the stochastic density fluctuations in the disk. 

Methods. We have performed 2D numerical simulations using a grid-based hydrodynamical code in which turbulence is modelled as stochastic 
i ■ forcing. We assume that the outermost planet is initially located just outside the position of the 3:2 mean motion resonance with the inner one 
and we study the dependance of the resonance stability with the amplitude of the stochastic forcing. 
~^ , Results. For systems of two equal-mass planets we find that in disk models with an effective viscous stress parameter a ~ 10~ 3 , damping 
^ effects due to type I migration can counteract the effects of diffusion of the resonant angles, in such a way that the 3:2 resonance can possibly 
remain stable over the disk lifetime. For systems of super-Earths with mass ratio q = m ; /m < 1/2, where m,(m ) is the mass of the innermost 
(outermost) planet, the 3:2 resonance is broken in turbulent disks with effective viscous stresses 2 x 10~ 4 < a < 1 x 10~ 3 but the planets become 
locked in stronger p + hp resonances, with p increasing as the value for a increases. For a > 2 x 10~ 3 , the evolution can eventually involve 
temporary capture in a 8:7 commensurability but no stable MMR is formed. 
^Tj ■ Conclusions. Our results suggest that for values of the viscous stress parameter typical to those generated by MHD turbulence, MMRs between 
two super-Earths are likely to be disrupted by stochastic density fluctuations. For lower levels of turbulence however, as is the case in presence 
of a dead-zone, resonant trapping can be maintained in systems with moderate values of the planet mass ratio. 
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1. Introduction PSR B1257+12, HD 69830, GJ 581 and HD 40307. For 

the sytems around main-sequence stars (HD 69830, GJ 581, 

To date, about 25 extrasolar planets with masses less than HD ^ obseryed period ^ between twQ adjacent 

10 and commonly referred to as super-Earths have been , , .. e e ... 

w J il low-mass planets are quite far from stnct commensurability. 

discovered (e.g. http://exoplanet.eu). Although two of tt , tU » • ,, , 

v 6 1 s—U. *L fe However, the planetary system that is orbiting the radio pulsar 

them, Corot-7b (Leger et al. 2009; Queloz et al. 2009) and GJ n( , D „,»,,,,» . , . ... ~ n ■ 

v 6 ' ^ ' PSR B 1257+ 12 exhibits two planets with masses 3.9 M B and 

1214b (Charbonneau et al. 2009) were detected via the transit A -, ■ -, .. IV , . . m , 

v ' 4.3 M© in a 3:2 mean motion resonance (Konacki & Wolszczan 

method, most of them were found by high-precision radial 2003) Papaloizou & Sz uszkiewicz (2005) showed that, for this 

velocity surveys. It is expected that the number of observed system ^ ^ existence of such a resonance can be understood by 

super-Earths will considerably increase in the near future with , , . , . , , . ... , 

F 1 a model in which two low-mass planets with mass ratio close 

the advent of space observatories Corot and Kepler. . .. , T • .. , nr , inm 

F F to unity undergo convergent type I migration (e.g. Ward 1997; 

Interestingly, Kepler team has recently announced the dis- rp , . , OAAO n . •, t U , , ■ , • 

6 -" F 1 Tanaka et al. 2002) while still embedded in a gaseous laminar 

covery of ~ 170 multi-planet systems candidates (Lissauer et disk ^ capture fa ^ resonance occurs More gen erall y , 

al. 2011), although these need to be confirmed by follow-up ^ amhors found ^ for more disparat£ mass ratios and 

programs. Previous to Kepler results, four multi-planet systems •■ ■ . . ... , .. e 

v fe provided that convergent migration occurs, the evolution of a 
containing at least two super-Earths had been detected around 
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system of two planets in the 1 - 4 M @ mass range is likely to 
result in the formation of high first-order commensurabilities 
p + hp with p>3. Studies aimed at examining the interaction 
of many embryos within protoplanetary disks also suggest 
that capture in resonance between adjacent cores through 
type I migration appears as a natural outcome of such a 
system (McNeil et al. 2005; Cresswell & Nelson 2006). This, 
combined with the fact that the majority of super-Earths are 
found in multiplanetary systems (Mayor et al. 2009), would 
suggest that systems of resonant super-Earths are common. 
The fact that most of the multiple systems of super-Earths 
observed so far do not exhibit mean motion resonances may 
be explained by a scenario in which strict commensurability 
is lost due to circularization through tidal interaction with the 
central star as the planets migrate inward and pass through the 
disk inner edge (Terquem & Papaloizou 2007). 

Moreover, it is expected that in presence of strong disk tur- 
bulence, effects arising from stochastic density flucuations will 
prevent super-Earths from staying in a resonant configuration. 
It is indeed now widely accepted that a source of anomalous 
viscosity due to turbulence is required to account for the 
estimated accretion rates for Class II T Tauri stars, which are 
typically ~ 10~ 8 Moyr 1 (Sicilia-Aguilar et al. 2004). The 
origin of turbulence is believed to be related to the magneto- 
rotational instability (MRI, Balbus & Hawley 1991) for which 
a number of studies (Hawley et al. 1996; Brandenburg et 
al. 1996) have shown that the non-linear outcome of this 
instability is MHD turbulence with an effective viscous stress 
parameter a ranging between ~ 5 x 10~ 3 and ~ 0.1, depending 
on the magnetic field amplitude and topology. 
So far, the effects of stochastic density fluctuations in the 
disk on the evolution of two-planet systems has received little 
attention. Rein & Papaloizou (2009) developed an analytical 
model and performed N-body simulations of two-planet 
systems subject to external stochastic forcing and showed that 
turbulence can produce systems in mean motion resonance 
with broken apsidal corotation, explaining thereby the resonant 
configuration of the HD 128311 system. Adams et al. (2008) 
examined the effets of turbulent torques on the survival of res- 
onances using a pendulum model with an additional stochastic 
forcing term. They found that mean motion resonances are 
generally disrupted by turbulence within disk lifetimes. 
Lecoanet et al. (2009) extended this latter work by considering 
disk-induced damping effets and planet-planet interactions. 
They found that systems with sufficiently large damping can 
maintain resonances and suggested that two-planet systems 
composed of a Jovian outer planet plus a smaller inner planet 
are likely to remain bound in resonance. 

In this paper we present the results of hydrodynamical 
simulations of systems composed of two planets in the 1 - 4 
M ffi mass range embedded in a protoplanetary disk in which 
turbulence is driven by stochastic forcing. Planets undergo 
convergent migration as a result of the underlying type I 
migration and we consider a scenario in which the initial 
separation between the planets is slightly larger than that 
corresponding to the 3:2 resonance. The aim of this work is to 



investigate whether or not resonant trapping can occur and be 
maintained in turbulent disks and how the stability of the 3:2 
resonance depends on the amplitude of the turbulence-induced 
density fluctuations. We find that for systems of equal-mass 
planets the 3:2 resonance can be maintained provided that 
the level of turbulence is relatively weak, corresponding to a 
value for the effective viscous stress parameter of a < 10~ 3 . 
In models with mass ratios q = m/m < 1/2 however, 
where m, (m ) is the mass of the inner (outer) planet, the 3:2 
resonance is disrupted in presence of weak turbulence but the 
planets can become eventually locked in higher first-order 
commensurabilities. For a level of turbulence corresponding 
to a ~ 5 X 10~ 3 however, MMRs are likely to be disrupted by 
stochastic density fluctuations. 

This paper is organized as follows. In Sect. 2, we describe 
the hydrodynamical model and the numerical setup. In Sect. 3, 
we use a simple model to estimate the critical level of turbu- 
lence above which the 3:2 resonance would be unstable. In 
Sect. 4 we present the results of our simulations. We finally 
summarize and draw our conclusions in Sect. 5. 



2. The hydrodynamical model 

2.1. Numerical method 

In this paper, we adopt a 2D disk model for which all the 
physical quantities are vertically averaged. We work in a 
non-rotating frame, and adopt cylindrical polar coordinates 
(R, <f>) with the origin located at the position of the central 
star. Indirect terms resulting from the fact that this frame is 
non-inertial are incorporated in the equations governing the 
disk evolution (Nelson et al. 2000). Simulations were per- 
formed with the FARGO and GENESIS numerical codes. Both 
codes employ an advection scheme based on the monotonic 
transport algorithm (Van Leer 1977) and include the FARGO 
algorithm (Masset 2000) to avoid timestep limitation due to 
the Keplerian orbital velocity at the inner edge of the grid. 
The evolution of each planetary orbit is computed using a 
fifth-order Runge-Kutta integrator (Press et al. 1992) and by 
calculating the torques exerted by the disk on each planet. 
We note that a softening parameter b = 0.6H, where H is the 
disk scale height, is employed when calculating the planet 
potentials. 

The computational units that we adopt are such that the 
unit of mass is the central mass M„, the unit of distance is the 
initial semi-major axis a, of the innermost planet and the unit 
of time is (GM./a 3 ) -1 ' 2 . In the simulations presented here, we 
use Nr = 256 radial grid cells uniformly distributed between 
Ri n = 0.4 and R ou , = 2.5 and N$ = 768 azimuthal grid cells. 
Wave-killing zones are employed for R < 0.5 and R > 2.1 in 
order to avoid wave reflections at the disk edges (de Val-Borro 
et al. 2006). 

In this work, turbulence is modelled by applying at each 
time-step a turbulent potential <b tur b to the disk (Laughlin & al. 
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2004, Bamteau & Lin 2010) and corresponding to the superpo- Model m, (M e ) m (M e ) I h 



sition of 50 wave-like modes. This reads: 
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with: Table 1. Parameters used in the simulations 



A* = &e 



cos(m k (p-(f> k -Cl k t k )sm(7Tf k /At k ). (2) 2.2. Initial conditions 



In Eq.|2] & is a dimensionless constant parameter randomly 
sorted with a Gaussian distribution of unit width. R k and <p k 
are, respectively, the radial and azimuthal initial coordinates 
of the mode with wavenumber m k , <r k = nR k IAm k is the radial 
extent of that mode and Q. k denotes the keplerian angular 
velocity at R — R k . R k and (f> k are both randomly sorted with 
a uniform distribution whereas m k is randomly sorted with 
a logarithmic distribution between m k = 1 and m k = 96. 
Following Ogihara et al. (2007), we set A* = if m k > 6 in 
order to save computing time. Each mode of wavenumber m k 
starts at time t = fo,* an d terminates when t k - t - fo,j > At k , 
where At k = 0.2nR k /m k c s , c s being the sound speed, denotes 
the lifetime of mode with wavenumber m k . Such a value for At k 
yields a turbulence with autocorrelation timescale t c ~ T or b, 
where T or t is the orbital period at R = 1 (Baruteau & Lin 
2010). 

In Eq. [1] y denotes the value of the turbulent forcing param- 
eter which controls the amplitude of the stochastic density 
perturbations. In the simulations presented here, we used 
four different values for y namely: y — 6 x 10~ 5 , 1.3 x 10~ 4 , 
1.9 x 10~ 4 , 3 x 10~ 4 . Given that y is related to the effective 
viscous stress parameter a and the disk aspect ratio h = H/R 
by the relation a = 1.4 x I0 2 (y/h) z (Baruteau & Lin 2010), the 
latter values for y correspond to a = 2 x 10~ 4 , 10~ 3 , 2 x 10~ 3 , 
5 x 10~ 3 respectively. Inviscid simulations with a — were 
also performed for comparison. 

In calculations with high values of y, viscous stresses aris- 
ing from turbulence can eventually lead to a significant change 
in the disk surface density profile over a few thousand orbits. 
This is also observed in three dimensional MHD simulations 
in which turbulence is generated by the MRI (Papaloizou & 
Nelson 2003). For lower values of y, such an effect also oc- 
curs but over a much longer timescale. In order to examine 
how this affects the results of the simulations, we have per- 
formed additional simulations in which the initial surface den- 
sity profile is restored on a characteristic timescale x m . We fol- 
low Nelson & Gressel (2010) and solve the following equation 
at each timestep: 

d~L _ X — 2,„,7 
dt t,„ 



(3) 



where E,„„ is the initial disk surface density and where r,„ was 
set to t,„ = 20 orbits. Such a value is shorter than the vis- 
cous timescale but longer than both the dynamical timescale at 
the outer edge of the disk and the lifetime of the mode with 
wavenumber m — 1 . The results of such simulations are dis- 
cussed in Sect. 14.1.41 



In this paper, we adopt a locally isothermal equation of 
state with a fixed temperature profile given by T — TqR~P 
where — 1 and where Tq is the temperature at R = 1 . This 
corresponds to a disk with constant aspect ratio h and for most 
of the simulations, we choose Tq so that h = 0.05. The initial 
surface density profile is chosen to be S,„,,(7?) = 'LoR~ cr with 
a = 0.5 and we have performed simulations with So = 2 x 10~ 4 
and So = 4 x 10~ 4 . Assuming that the radius R = 1 in the 
computational domain correponds to 5.2 AU, such values 
for So correspond to disks containing 0.02 M* and 0.04 M* 
respectively of gas material interior to 40 AU. No kinematic 
viscosity is employed in all the runs presented here. 

The inner and outer planets initially evolve on circular or- 
bits at a, = 1 and a = 1.33 respectively, which corresponds 
to a configuration for which the outermost planet is initially lo- 
cated just outside the 3:2 MMR with the inner one. For most 
models, we focus on equal-mass planets with m, = m < 3.3 
M e , where m, (m ) is the mass of innermost (outermost) planet. 
However, we have also considered one case in which the planet 
mass ratio q = mijm v&q- 1/2. The parameters for all models 
we conducted are summarized in Table Q] Given that the type I 
migration timescale T mlgjP of a planet with mass m p , semimajor 
axis a p and on a circular orbit with angular frequency Q. p can 
be estimated in the locally isothermal limit by (Paardekooper 
et al. 2010): 



1 mtg,p 



(1.6+J0 + O.7O-) 



— — ^-t^q: 1 , 

m p T.(a p )rj 



(4) 



we expect that equal low-mass planets embedded in our disk 
model will undergo convergent migration and become eventu- 
ally trapped in the 3:2 resonance. For a larger initial separation 
between the two planets, capture in 2: 1 resonance may also oc- 
cur. However, test simulations have shown that unless the disk 
mass is very low, differential migration is not slow enough for 
the planets to become trapped in that resonance. This justifies 
our assumption that the planets are initially located just outside 
the 3:2 resonance. We also comment that equal-mass planets 
migrating in the type I regime will undergo convergent migra- 
tion provided that cr < 3/2 whereas cr > 3/2 will lead to diver- 
gent migration. 

3. Theoretical expectations 

In this section, we consider two low-mass planets embedded in 
a turbulent disk and locked in a p + 1 :p mean motion resonance, 
and we derive the critical amplitude of the turbulent forcing 
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Fig. 1. Time evolution of the period ratio resulting from N-body 
runs for model Gl and for six different realizations with y = 
1.3 x l(T 4 ,y = 1.9 x 1(T 4 , andy = 3 x 1(T 4 . 



y c below which the resonance would remain stable. Following 
Adams et al. (2008) and Rein & Papaloizou (2009), we assume 
that only the outermost planet experiences the torques arising 
from the disk. We also assume that the planets have near equal 
mass, in order to avoid the chaotic regime which comes into 
play for disparate masses (Papaloizou & Szuszkiewicz 2005). 
In the limit of large damping rate for the resonance and neglect- 
ing effects from planet-planet interaction, the asymptotic value 
P of the probability for the resonance to be maintained is given 
by (Lecoanet et al. 2009): 



p = 4 



1 



1/2 



(5) 



where is the diffusion coefficient associated with the res- 
onant angle diffusion and is the damping timescale for the 
resonance angle. This equation is valid at late times t » a> Q l 
where a>o is the libration frequency of the resonant angles, as 
long as t » D^ 1 and D^ 1 » tj. From the previous equation, we 
can estimate the maximum value of the diffusion coefficient for 
the sytem to remain bound in resonance with probability P = 1 . 
This reads: 



16 
= — t , 



(6) 
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Fig. 2. Upper panel: time evolution of the resonant angles 0i = 
3A -2Ai-a>i (black) and cp2 = 3A -2Ai- a> (red) resulting from 
N-body runs for model Gl and for two different realizations 
with y = 1 .3 x 10~ 4 . Middle panel: same but for y = 1 .9 x 10~ 4 . 
Lower panel: same but for y = 3 x 10" 4 . 

As shown in Adams et al. (2008), can be related to the diffu- 
sion coefficient Dh,o associated with the diffusion of the outer 
planet's angular momentum as: 



Dh,o 



For moderate eccentricities, a>o is given by: 



with C = q Q.iaf d (a), 



(7) 



(8) 



where e-, is the eccentricity of the inner planet, go = mo/M* 
and Q,- is the angular frequency of the innermost planet. In 
the previous equation, a = a,i/a , (j2,ji) are integers which 
depend on the resonance being considered and fd(a) results 
from the expansion of the disturbing function. In the case of the 

3:2 resonance, we have 72 = -2, 74 = -1 and afd{a) 1.54 

(Murray & Dermott 1999). 

In Eq.|7l Dh,o can be expressed in terms of both the correla- 
tion timescale t c associated with the stochastic torques exerted 
on the outer planet and the standard deviation of the turbulent 
torque distribution o~ t as: 



D 



H,o 



(9) 



As discussed in Baruteau & Lin (2010), cr, takes the fol- 
lowing form when applied to the outermost planet: 
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cr, = CL q ya 4 Q. 2 , 

where q = ra /M+, E is the value of the surface density 
at the position of the outer planet, D. is the angular frequency 
of this planet and C is a constant. For a simulation using the 
same disk parameters as for model Gl and with y — 6 x 10~ 5 , 
we find C ~ 1.6 X 10 2 , which is close to the value found by 
Baruteau & Lin (2010). Combining Eqs. 171 19"1 and [TOl gives an 
expression for the diffusion coefficient associated with the dif- 
fusion of the resonant angle D®, in terms of the value for the 
turbulent forcing y. This reads: 



D„ = 



2C 2 q yn 3 



(11) 



with qd = TiLocrJMi,. Setting 6a> = ojq/Q. , we can rewrite the 
previous expression as: 

2c 2 q yn 



D,„ 



(12) 



9tt6u 2 

We notice that in the case where p — 1, 6a> is comparable to 
the dimensionless width of the libration zone. Using the previ- 
ous equation together with Eq. [6] we find that the critical value 
for the turbulent forcing above which the 3:2 resonance is dis- 
rupted is given by: 

7,- 5.3X10- 2 — (r d Q o r 1/2 . 
qd 

In absence of turbulent forcing, we expect the amplitude of 

— 1/2 

the resonant angles to scale as £2 G ' (Peale 1976). According 
to Rein & Papaloizou (2009), this implies that the damping 
timescale of the libration amplitude is twice the migration 
timescale T mlf of the whole system, namely that composed of 
the two planets locked in resonance and migrating inward to- 
gether. In that case, the previous equation becomes: 



5u> 



3.7 x 10- 2 — (r„„„Q ) 



■1/2 



qd 



In the limit where m* ~ 0, we would have r, 



where T m i gfi is the migration timescale of the outer planet which 
is given by Eq.|4]with h - 0.05, cr = 0.5 and (3 — 1. In that case 
the expression for y c becomes for our disk model: 

y c ~ 3.5 x \Q- 2 8(oq]j 2 q d XI2 hr l . 

In order to check the validity of the previous expression for 
7c, we have performed a few N-body runs using a fifth-order 
Runge Kutta method. In these calculations, the forces arising 
from type I migration are not determined self-consistently but 
modelled using prescriptions for both the migration rate r ap and 
eccentricity damping rate t C(j of the planets. For r a , we used: 

( l + (e P /h) 5 \ 

T " = [i-ie p /hy* )*•*•>' (16) 

where e p is the planet eccentricity and where T m i glP is given by 
Eq. |4] and where the numerical factor accounts for the modifi- 
cation of the migration rate at large eccentricities (Papaloizou 
& Larwood 2000). For T e> we used (Tanaka & Ward 2004): 



Te p = -^(l+0.25(e p /hf)T„ 



In the last equation, K ~ 1.7 is a constant which was cho- 
(10) sen in such a way that the eccentricity damping rate obtained 
in N-boby runs gives good agreement with that resulting from 
hydrodynamical simulations. Following Rein et al. (2010), we 
model effects of turbulence as an uncorrelated noise by per- 
turbing at each time step At the velocity components v; iP of 
each planet by Av,-. p = V2DAf£ where £ is a random variable 
with gaussian distribution of unit width. D is the diffusion coef- 
ficient which should vary as the planets migrate but which was 
fixed here to a constant value of D = crfrc/a 2 . 
In Fig. Q] we show the time evolution of the orbital period ra- 
tio for N-body simulations with parameters corresponding to 
model Gl and for six different realizations with y = 1.3 x 
10~ 4 , 1.9 x 10~ 4 ,3 x 10" 4 . For this model, we estimate 6u ~ 
1.9 x 10~ 3 (see Sect. H~L4l which leads to y c ~ 2.5 x 10~ 4 
using Eq. Q3] From Fig. [1] it appears that capture in the 3:2 
MMR occurs for most of the realizations with y < 1.9 x 10~ 4 . 
For two specific realizations of each value for y we considered, 
the time evolution of the resonant angles (p\ - 3A - 2A, - co, 
and <p2 - 3A„ - 2/1, - a> associated with the 3:2 resonance, 
where Ai (A ) and w,- (w ) are respectively the mean longitude 
and longitude of pericentre of the innermost (outermost) planet, 
is displayed in Fig. [2] Although the angles can eventually cir- 
culate for short periods of time, it is clear that the 3:2 MMR 
remains stable on average for y < 1 .9 x 10~ 4 . For y — 3 x 10~ 4 , 
(13) we find that the planets pass through the 3:2 resonance in two 
of the six realizations performed while the four other can even- 
tually involve temporarily capture in the 3:2 MMR. In these 
cases however, the lifetime of the resonance does not exceed a 
few hundred orbits, as can be seen in the lower panel of Fig. [2] 
which displays for y — 3 x 10~ 4 the time evolution of tf>\ and 
<p2 for two realizations in which the period ratio remains close 
to that corresponding to the 3:2 MMR. Therefore, the results 
of these N-body calculations suggest that the 3:2 MMR is only 
marginally stable for such a value of y, which is consistent with 
the aforementioned analytical estimation of y c ~ 2.5 x 10~ 4 . 

4. Results of hydrodynamical simulations 

For equal mass planets (q — 1), results of hydrodynamical sim- 
ulations suggest that capture in 3:2 resonance can occur in tur- 
n 5) bulent disks for which the level of turbulence is relatively weak. 
For systems with q < 1/2 however, it appears that trapping in 
the 3:2 resonance is maintained only provided that the disk is 
close to being inviscid. 

4.1. Models with q - 1 

For inviscid simulations with equal low-mass planets, the abil- 
ity for the two planets to become trapped in the 3:2 resonance 
depends mainly on the planets' relative migration rate which 
scales as hr 2 . For model G3 (h = 0.04), we find that capture 
in 3:2 resonance does not occur in that case due to the relative 
migration timescale being shorter than the libration period cor- 
responding to that resonance. For other models with h = 0.05 
however, it appears that the system can enter in a 3:2 commen- 
surability which remains stable for the duration of the simula- 
(17) tion, which generally covers ~ 10 4 orbits at R — 1. This occurs 



(14) 
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Fig. 3. Upper left (first) panel: time evolution of planet semi-major axes for model Gl and for the different values of y we 
considered namely for y — (black line), y = 6 x 1CT 5 (red line), y = 1 .3 x 1CT 4 (blue), y = 1 .9 x 1CT 4 (green) and y = 3 x 1CT 4 
(orange). Upper right (second) panel: time evolution of the inner planet eccentricity. Third panel: time evolution of the outer 
planet eccentricity. Fourth panel: time evolution of the period ratio. Simulations were performed with GENESIS. 



not only for laminar disks, but also for turbulent disks provided 
that the value for the turbulent forcing is not too large. For ex- 
ample, we find that the 3:2 commensurability is maintained in 
most of the turbulent runs with y < 1 .3 x 10~ 4 in models Gl and 
G2 whereas in model G4 this occurs provided that y < 6x 10 5 . 
Below we describe in more details the results of the simulations 
with q = 1 and we use model Gl to illustrate how the evolution 
depends on the value for the forcing parameter y. 

4.1 .1 . Orbital evolution 

The time evolution of the planets' semi-major axes, eccentrici- 
ties and period ratio corresponding to model Gl and for one re- 
alization of the different values of y we considered are depicted 
in Fig. [3] In each case, the period ratio is observed to initially 
decrease, suggesting that early evolution involves convergent 
migration of the two planets. Not surprisingly, a tendancy for 
the planets to undergo a monotonic inward migration is ob- 
served at the beginning of the simulations with the lowest val- 
ues of y whereas these are more influenced by stochastic forc- 
ing for y > 1 .9 x 10 4 . This is due to the fact that the amplitude 
of the turbulent density fluctuations is typically stronger than 
that of the planet's wake for simulations with y > 1.9 x 1CT 4 , 
as shown in Fig. [4] which displays snapshots of the perturbed 
surface density of the disk for different values of y. 
The semimajor axis evolution also reveals a clear tendancy 
for lower migration rates with increasing y, which is an effect 



arising from the desaturation of the horseshoe drag by turbu- 
lence (Baruteau & Lin 2010). As y increases, the disk torques 
are indeed expected to increase from the differential Lindblad 
torque, obtained for y — 0, up to the fully unsaturated torque, 
which is maintained for a ~ 0.16(m,/M») 3 / 2 /i~ 4 (Baruteau & 
Lin 2010). For h = 0.05, such a value for a corresponds to 
y ~ 1.2 x 10~ 4 . For higher values of y however, we expect the 
torques to slightly decrease with increasing y due to a cut-off of 
the horseshoe drag arising when the diffusion timescale across 
the horseshoe region is smaller than the horseshoe U-turn time 
(Baruteau & Lin 2010). 

This can be confirmed by inspecting the evolution of the 
running-time averaged torques exerted on both planets and 
which are presented in Fig. [5] Up to a time of approximately 
~ 10 3 orbits and for y < 1.3 x 10~ 4 , the torques are observed to 
increase with increasing y, while they decrease for higher val- 
ues of y. This is in very good agreement with the expectation 
that the fully unsaturated torque is reached for y — 1.2 x 10~ 4 . 
At later times however, the torques obtained in simulations with 
y > 1.9 x 10~ 4 can eventually exceed those computed in runs 
with y < 1.3 x 10~ 4 . We suggest that this is related to the fact 
that the disk surface density tends to be significantly modified 
at the planet positions at high turbulence level. This is illus- 
trated in Fig. [6] which shows the disk surface density profile 
at t = 2000 orbits for the different values of y we consid- 
ered. Here the inner and outer planets are located at a, ~ 0.98 
and a Q ~ 1.25 respectively. It is interesting to note that for 
y — 3 x 10~ 4 , the outer planet tends to evolve in a region of 
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Fig. 4. This figure shows, for model Gl, snapshots of the perturbed surface density of the disk for y — (first panel), y = 6 x 10 5 
(second panel), y = 1.3 x 10~ 4 (third panel) and y = 1.9 x 10~ 4 (fourth panel). 
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Fig. 5. Upper panel: time evolution of the running-time av- 
eraged torques exerted on the inner planet for model Gl and 
for the different values of y we considered namely for y = 
(black line), y = 6 x 1(T 5 (red line), y = 1.3 x 1(T 4 (blue), 
y = 1 .9 x 10~ 4 (green) and y = 3 x 10~ 4 (orange). Lower panel: 
same but for the outer planet. 



Fig. 6. Disk surface density profile at t = 2000 orbits for model 
Gl and for the different values of y we considered namely for 
y = (black line), y = 6 x 10~ 5 (red line), y = 1.3 x 10~ 4 
(blue), y = 1.9 x 10~ 4 (green) and y = 3 x 10~ 4 (orange). 



positive surface density gradient where the corotation torque is 
positive, in such a way that the torque exerted on that planet 
can become higher than that exerted on the inner one. This ef- 
fect is responsible for the increase of period ratio observed at 
late times in simulations with y > 1.9 x 10~ 4 . 
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Fig. 7. Upper panel: time evolution of the theoretical change of 
the inner planet eccentricity de,jdt given by Eq.[18]for model 
Gl and for the different values of y we considered namely for 
y = (black line), y = 6 x 1(T 5 (red line), y = 1 .3 x 1(T 4 (blue), 
y = 1.9 x 1CT 4 (green) and y — 3 x 1CT 4 (orange). Lower panel: 
same but for y — 6 x 10~ 5 and models Gl (red), G2 (blue) and 
G4 (green). 



4.1.2. Eccentricity evolution 

For model Gl, examination of the early planets' eccentricities 
evolution displayed in Fig.[3]shows a clear tendancy of higher 
eccentricities with increasing the value for y, which is in agree- 
ment with the expectation that turbulence is a source of ec- 
centricity driving (Nelson 2005). This can be confirmed by in- 
specting the theoretical change of the inner planet eccentricity 
dei/dt which can be computed using the following expression 
(Bums 1976): 



de, 
dt 



2e? 



E 

- +2— 

E H 



(18) 



where E 
inner planet, H 



-G(Mi, + mi) 1 2a i is the specific energy of the 

= ^G{M+ + m;)fl;(l - e?) its specific angular 

momentum, H the torque exerted by the disk and E the power 
of the force exerted by the disk on the planet. The early time 
evolution of de-,jdt is displayed, for this model and for the dif- 
ferent values of y we considered in the upper panel of Fig. Q It 
clearly demonstrates that, compared with the laminar run, the 
theoretical rate of change of e, is higher in turbulent runs and 
that it increases with increasing the value for y. 
Also, inspection of the lower panel of Fig. |7]reveals that, com- 
pared with model Gl in which m, = m = 3.3 M ffi and 
So = 2 x 10~ 4 , the disk induced eccentricity damping is weaker 
in model G4 for which m, = m = 1.6 M e . This is due to the 
damping of eccentricity at coorbital Lindblad resonances scal- 
ing linearly with planet mass (Nelson 2005). Given that this 
damping also scales with disk mass, this explains why the disk 
induced eccentricity damping is stronger in model G2 in which 
S = 4 x 10" 4 . 
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Fig. 8. Upper panel: time evolution of the period ratio for 
model Gl and for four different realizations with y = 1 .3x 10 -4 . 
Lower panel: same but for y = 1.9 x 10~ 4 . Simulations were 
performed with GENESIS. 



4.1.3. Time evolution of the resonant angles 

As mentioned above, we find that for model Gl a stable 3:2 
commensurability generally forms for y < 1.3 x 10~ 4 , which 
can be confirmed by inspecting the upper panel of Fig. [8] 
which displays the time evolution of the period ratio for four 
different realizations with this value of y. However, as one can 
see in the lower panel of Fig. [8] such a resonance is observed 
to break at late times for all the realizations performed with 
y > 1.9 x 10 -4 . Comparing Fig. [8] and the two upper panels 
in Fig. [1] we see that the results of these hydrodynamical 
simulations are at first sight in good agreement with those of 
N-body runs. Moreover, it is interesting to note that in some 
runs, planets which leave the 3:2 resonance can eventually 
pass through that resonance again at later times or others 
commensurabilities like the 2:1 resonance. 
A similar outcome is observed in model G2 which had a value 
of Eo = 4 X 10~ 4 . This arises because both the turbulent torque 
and the damping rate of the libration amplitude scale linearly 
with So. For model G4 however, which had E u = 2 X 10 4 and 
nii - tn = 1.6 Me, the 3:2 resonance is disrupted in the run 
with y - 1.3 x 10~ 4 due to the damping rate being smaller 
for this model. For a single realization for each value of y, the 
evolution of the period ratio for models G2 and G4 is displayed 
inFig.g] 

In the two upper panels of Fig. [I0]we show, for model Gl 
and for two realizations with y - 6 x 10~ 5 and y = 1.3 x 10~ 4 , 
the time evolution of the resonant angles associated with the 3:2 



A. Pierens et al.: On the dynamics of resonant super-Earths in disks with turbulence driven by stochastic forcing 



9 



1.60 




6| r i"tij5fljl|iin<i<f , »" 



1000 2000 3000 4000 5000 6000 7000 

Time (orbits) 




6.0-1 J 8.0-1 J 
Time (orbits) 



1.0-10" 1.2-10" 1.4-10" 



Fig. 9. Upper panel: time evolution of the period ratio for 
model G2 and for the different values of y we considered 
namely for y = (black line), y = 6xl(T 5 (red), y = 1.3xl0~ 4 
(blue), y = 1.9 x 10~ 4 (green) and y = 3 x 10~ 4 (orange). Lower 
panel: same but for model G4. Simulations were performed 
with GENESIS. 

resonance. For y — 6 x 1CT 5 , the 3:2 resonance is established at 
t ~ 1800 orbits while it forms at t ~ 2000 orbits for the calcu- 
lation with y = 1.3 x 10~ 4 . This is consistent with the fact that, 
for moderate values of y, migration rates tend to decrease with 
increasing y. As illustrated in the second and third panels of 
Fig. [3] resonant capture makes the eccentricities of both plan- 
ets grow rapidly before they saturate at values of <?,■ ~ e a ~ 0.01 
in the run with y — 6 x 10~ 5 . As discussed in Sect. 14.1.21 these 
tend to be larger in the case where y — 1 .3 x 10~ 4 , with the ec- 
centricities reaching peak values of e, ~ 0.02 and e a ~ 0.015. 
Not surprisingly, there is a clear trend for the amplitude of 
the resonant angles to increase with increasing the value for 
y in model Gl. For the run with y — 6 x 10~ 5 , the angles 
slightly spread until t ~ 5 x 10 3 orbits before their ampli- 
tude continuously decrease with time. This indicates that over 
long timescales damping of the resonant angles through mi- 
gration tends to overcome diffusion effects. In the case where 
y = 1.3 x 10~ 4 however, periods of cyclic variations of the 
resonant angles can be seen with the angles librating with high 
amplitude before being subsequently damped. Given that in ab- 
sence of turbulent forcing, the libration amplitude should de- 

— 1/2 

crease as Q, (Peale 1976), we would expect the 3:2 reso- 
nance to be maintained, for y < 1.3 x 10~ 4 , over timescales 
much longer than those covered by the simulations. 
For two realizations with y = 6 x 10~ 5 and y = 1 .3 x 10~ 4 , the 
evolution of both 4>\ and 4>2 for models G2 and G4 is diplayed 
in the middle and lower panels of Fig. [I0]respectively. In com- 
parison with model Gl, the resonant angles librate with slightly 
higher amplitudes in model G2 and can eventually start oscil- 



5. 

2 
1 



6 - 
5 

* 3 
2- 

1 



W0. 



ii,i;f : 



• *B:WVI , , , 

2000 4000 6000 8000 10000 2000 4000 6000 8000 10000 
Time (orbits) Time (orbits) 




6 r 
5 - 
4 

3 y- 

2 - 
1 % 



:•: < f i 



|i 1 1 



\ MM 



1400 2800 4200 5600 7000 
Time (orbits) 



1400 2800 4200 5600 7000 
Time (orbits) 




6 
5 

4- 
3 
2- 
1 



5.0-1 3 1.0-10" 1.5-10" 

Time (orbits) 



5.0-1 3 1.0-10" 1.5-10" 
Time (orbits) 



Fig. 10. Upper panel: time evolution of the resonant angles 
<pi = 3A - 2Ai - iOj (black) and 02 = 3/t„ - 2Aj - w„ (red) 
for model Gl with y = 6 x 10~ 5 (left) and y = 1.3 x 10 4 
(right). Middle panel: same but for model G2. Lower panel: 
same but for model G4. 



lations between periods of circulation and libration in the run 
with y — 1.3 x 10~ 4 . Again this arises because, compared with 
model Gl, the turbulent density fluctuations are stronger in this 
model. In model G4 however, the 3:2 resonance is maintained 
for only ~ 3 x 10 3 orbits in the case where y — 1 .3 x 10 4 , which 
indicates that for this model the damping rate is too weak for 
this resonance to remain stable. 

4.1.4. Comparison with analytics 

We now examine how the results of the simulations described 
above compare with the expectations discussed in Sect. [3] For 
model Gl, we can estimate the libration frequency oiq using 
Eq.[8]in conjunction with the results from this model shown in 
Fig. [3] Adopting the inviscid simulation as a fiducial case, we 
have a, = 0.9, a = 111 and e, = 0.01 at t ~ 4000 orbits, 
which leads to 6a> = u>o/Q. ~ 1.9 x 10~ 3 . Moreover, the migra- 
tion timescale of the system was estimated to T m ; g ~ 3.3 x 10 4 
orbits from the results of this simulation. Using Eq. H4lwhich 
gives the value for y c as a function of T mlg , we can then pro- 
vide an analytical estimate y™° of the critical amplitude for 
the turbulent forcing above which the 3:2 resonance should be 
disrupted. For this model, this critical value is estimated to be 
yf ~ 1.9 x 10~ 4 while Eq. [B] predicts y a c " a ~ 2.5 x 10~ 4 . 
Returning to Fig. [8] we see that the results of the simulations 
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Fig. 11. Upper panel: time evolution of the period ratio for 
model Gl and for four runs with y — 6 x 10~ 5 and in which 
Eq. [3] is solved at each timestep . Lower panel: same but for 
y — 1.3 x 10~ 4 . Simulations were performed with FARGO. 



in which a roughly constant surface density profile is main- 
tained (see Sect. 12. U produced slightly different results since 
we find 6 x 1CT 5 < y c < 1.3 x 1CT 4 in these cases. The little 
difference exhibited by our two codes is apparently due to the 
fact that turbulence induces changes in the surface density pro- 
file which are slightly different. In FARGO, the disk density at 
the position of the inner planet is slightly higher compared with 
GENESIS while it is slightly lower at the position of the outer 
planet. 

For calculations in which Eq.[3]is solved at each timestep, the 
time evolution of the period ratio for four realizations with 
y = 6 x 1(T 5 and y = 1.3 x 1(T 4 is displayed in Fig.Ql] In that 
case, all the realizations performed with y = 6x10 5 resulted in 
the formation of the 3:2 resonance whereas for y — 1.3 x 1CT 4 , 
two of the four realizations resulted in capture in that resonance 
by the end of the run. For y - 1.3 x 1CT 4 , the time evolu- 
tion of the resonant angles associated with the 3:2 resonance 
is displayed in Fig. [12] Compared with previous runs in which 
the surface density profile was altered by turbulence, we see 
that capture in resonance tends to occur later in runs where a 
roughly constant surface density profile is maintained. This oc- 
curs because the disk density at the position of the inner planet 
tends to be higher in runs where the initial disk surface density 
profile is restored, leading to a slower differential migration be- 
tween the two planets. 

For other models, repeating the previously decribed procedure 
leads to analytical estimates of y™" = 1.2 x 10~ 4 for model G2 
and y" na ~ 9.3 x 1(T 5 for model G4. Given that the simula- 
tions performed with GENESIS indicate that 1.3 x 10~ 4 < y c < 
1 .9 x 1 - 4 and 6 x 1 (T 5 < y c < 1 .3 x 1 (T 4 for models G2 and G4 
respectively, we see that again the previous analytical estimates 
compare reasonably well with the results of our simulations. 



performed with GENESIS suggest that 1.3 x 1(T 4 < y c < 
1.9 x 1CT 4 for this model, which is clearly in broad agreement 
with the previous analytical estimate. We note however that 
both simulations performed with FARGO and additional runs 
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Fig. 12. Time evolution, for model Gl and for four different 
realizations, of the resonant angles 0i = 3/l e - 2/1, - w, (black) 
and <p2 = 3A a - 2Ai - <±> (red) for y = 1 .3 x 1CT 4 and in the case 
where Eq.|3]is solved at each timestep. 



4.2. Model with q= 1/2 

For systems with q — 1/2, the results of the simulations in- 
dicate that the 3:2 resonance can be maintained only in cases 
where the disk is close to being inviscid. Fig.[T3lshows the re- 
sults for model G5 in which m, = 1.6 M® and m - 3.3 M ffi 
and for a single realization of the different values of y we con- 
sidered. Moving from left to right and from top to bottom the 
panels display the time evolution of the planets' semi-major 
axes, eccentricity of the inner planet, eccentricity of the outer 
one, and period ratio. The evolution of semi-major axes shows 
strong similarities with that of model Gl, with the migration 
rates of the planets observed to decrease with increasing the 
value for y, as discussed in Sect. 14.1.11 

Because of stochastic density fluctuations, the planets eccen- 
tricities are highly variables quantities and a clear trend of 
higher eccentricities for higher values of y is again observed 
at the beginning of the simulations. For y — 0, the time evolu- 
tion of the period ratio shows that trapping in the 3:2 resonance 
occurs at t ~ 10 3 orbits. This resonant interaction causes ec- 
centricity growth of both planets with the eccentricities of the 
inner and outer planets reaching peak values of e, ~ 0.035 and 
e a ~ 0.02 respectively, although convergence is not fully es- 
tablished at the end of the simulation. Comparing Figs [3] and 
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Fig. 13. Upper left (first) panel: time evolution of planet semi-major axes for model G5 and for the different values of y we 
considered namely for y — (black line), y — 6 x 1CT 5 (red line), y = 1.3 x 1CT 4 (blue), y = 1.9 x 10~ 4 (green) and y = 3 x 1CT 4 
(orange). Upper right (second) panel: time evolution of the inner planet ecentricity. Third panel: time evolution of the outer 
planet eccentricity. Fourth panel: time evolution of the period ratio (a /a,) 15 . Simulations were performed with GENESIS. 



[T3l we see that the period ratio oscillates with greater ampli- 
tude in model G5, indicating thereby that resonant locking is 
weaker. This arises because, compared with model Gl, the rel- 
ative migration rate is higher for that model. Here, it is worth- 
while to notice that over timescales longer that those covered 
by the simulations, it is not clear whether or not the planets will 
remain bound in the 3:2 resonance since for disparate planet 
masses as is the case for model G5, we expect the dynamics 
of the system to be close to the chaotic regime even for y = 
(Papaloizou & Szuszkiewicz 2005). 

For turbulent runs, we find that the planets become temporarily 
trapped in the 3:2 resonance but in each case, the final out- 
come appears to be disruption of that resonance. Not surpris- 
ingly, the survival time of the 3:2 resonance tends to increase 
with decreasing the value for y. For these realizations, the life- 
times of the resonance are estimated to be ~ 1500 orbits for 
y = 6 x 10~ 5 , ~ 2 x 10 3 orbits fory = 1.3 x 10~ 4 , ~ 10 3 orbits 
for y = 1 .9 x 10~ 4 . The slightly longer lifetime of the resonance 
obtained in the run with y = 1.3 x 10~ 4 arises because of the 
stochastic nature of the problem. 

In Fig. [14] is displayed for two runs with y = 6 x 10 and 
y = 1.3 x 10~ 4 the evolution of the resonant angles associated 
with the 3:2 resonance. In comparison with models in which 
q = 1, there is a clear tendancy for these to librate with higher 
amplitude. Therefore, for models with disparate planet masses, 
the origin of the disruption of the 3:2 resonance in turbulent 
runs is likely to be related to the combined effect of diffusion 
of the resonant angles plus high libration amplitudes due to the 



resonance being weaker. 

For y > 6 x 10~ 5 , we can not rule out the possibility that 
the planets would become locked in stronger p + l:p reso- 
nances with p > 3. To investigate this issue in more details, 
we have performed an additional set of simulations in which, 
for each value of y we considered, the outer planet was ini- 
tially located just outside the 4:3 resonance with the inner one. 
In cases where the 4:3 resonance was found to be unstable, we 
performed an additional run but with an initial separation be- 
tween the two planets slightly larger than that corresponding 
to the 5:4 resonance. If the 5:4 resonance is not maintained, 
we repeat the procedure described above until a stable p + 1 :p 
commensurability is eventually formed. Because performing 
several realizations for each value of y would require a large 
suite of simulations, we have considered here only one single 
realization for each value of y. 

Fig. Q3] illustrates how the established resonance depends on 
the value for the turbulent forcing. Not surprisingly, a clear 
trend for forming stronger p + l:p resonances with increas- 
ing y is observed. For y = 6 x 10~ 5 , the system enters in a 
stable 4:3 resonance while for y = 1.3 x 10~ 4 , the planets 
become rather locked in the 5:4 resonance. For the runs with 
y = 1.9 x 10~ 4 and y = 3 x 10~ 4 however, the planets become 
temporarily trapped in the 8:7 resonance but in each case this 
commensurability is subsequently lost with the planets under- 
going divergent migration. In that case, it is interesting to note 
that the system is close to the stability limit since for planets in 
the Earth mass range as is the case here, we expect resonance 
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overlap to occur for p > 8 (Papaloizou & Szuszkiewicz 2005). 
Therefore, we can reasonably suggest that for such values of y, 
super-Earths with mass ratio q = mj/m < 1/2 may not be able 
to become trapped in a stable mean motion resonance and may 
eventually suffer close encouters. 

5. Discussion and conclusion 

In this paper we have presented the results of hydrodynamic 
simulations aimed at studying the evolution of a system com- 
posed of two planets in the Earth mass range and embedded 
in a turbulent protoplanetary disk. We employed the turbulence 
model of Laughlin et al. (2004) and modified by Baruteau & 
Lin (2010) in which a turbulent potential corresponding to the 
superposition of multiple wave-like modes is applied to the 
disk. We focused on a scenario in which the outermost planet is 
initially located just outside the 3:2 resonance and investigated 
how the evolution depends on both the planet mass ratio q and 
the value for the turbulent forcing parameter y. 
The results of the simulations indicate that for systems with 
equal mass planets, a 3:2 resonance can be maintained in pres- 
ence of weak turbulence. For instance, in the case of two plan- 
ets with equal mass 3.3 M ffi , we find that the 3:2 resonance is 
stable in runs with y < 1.9 x 10~ 4 , which corresponds to values 
for the effective viscous stress parameter of a < 2 x 10~ 3 . Such 
a value was found to compare fairly well with that resulting 
from both analytical estimations and preliminary N-body runs. 
For systems with planet mass ratios q < 1/2 however, it ap- 
pears that a 3:2 resonance can remain stable only provided that 
the disk is close to being inviscid. In turbulent disks however, 
the outcome depends strongly on the value for y: 

i) For 6 x 10~ 5 < y < 1 .3 x 10~ 4 (equivalent to 2 x 10~ 4 < a < 
10~ 3 ) , the planets tend to become locked in stronger p + l:p 
resonances, with p increasing as the value for y increases. 

ii) In the case where y > 1 .9x 10~ 4 (equivalent to a > 2x 10 -3 ), 
we find that the planets can become temporarily trapped in a 8:7 
commensurability, but this resonance is disrupted at later times 
and no stable resonance is formed. 

Given that the volume averaged stress parameter deduced from 
MHD simulations is typically a ~ 5 x 10~ 3 (Papaloizou & 
Nelson 2003; Nelson 2005), these results suggest that mean 
motion resonances between planets in the Earth mass range 
are likely to be disrupted in the active zones of protoplanetary 
disks. For relatively low levels of turbulence however, as is the 
case for a dead-zone (Gammie 1996), a resonance can be main- 
tained for moderate values of the planet mass ratio. 
Such a scenario is broadly consistent with the preliminary anal- 
ysis of ~ 170 multi-planetary systems candidates recently de- 
tected by Kepler (Lissauer et al. 201 1) and which suggests that 
only a few of the observed adjacent pairings are either in or near 
a MMR. However, examination of the slope of the cumulative 
distribution of period ratios (Fig. 7 of Lissauer et al. 201 1 ) also 
reveals an excess of planets with period ratios corresponding to 
the 2:1 or 3:2 commensurabilities. In that case, it appears that 
the neighboring planet candidates have masses within 20% of 
each other. This clearly supports our findings that in disks with 
moderate levels of turbulence, MMRS are stable provided the 
mass ratio between the neighboring planets is close to unity. 
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Fig. 14. Time evolution of the resonant angles <p[ = 3A - 2 A, - 
u>i (black) and <p2 = 3/l - 2/1,- - a> (red) for model Gl with 
y = 6 x 10~ 5 (left panel) and y = 1.3 x 10~ 4 (right panel). 
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Fig. 15. Period ratio between the two planets for model G5 with 
y = (black line), y = 6 x 10~ 5 (red), y = 1.3 x 10~ 4 (blue), 
y = 1.9 x 10~ 4 (green) and y = 3 x 10~ 4 (orange). Simulations 
were performed with GENESIS. 

Since turbulence has a significant impact on the capture of two 
planets in the Earth mass range, it will be of interest to exam- 
ine this issue using three-dimensional MHD simulations, which 
eventually include the presence of a dead-zone. We will address 
this issue in a future paper. 
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